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ABSTRACT 

Gravitational waves are predicted by Einstein's theory of general relativity as well as 
i-Q ■ other theories of gravity. The rotational stability of the fastest pulsars means that 

Mh' timing of an array of these objects can be used to detect and investigate gravitational 

Q , waves. Simultaneously, however, pulsar timing is used to estimate spin period, period 

$— ( ' derivative, astrometric, and binary parameters. Here we calculate the effects that a 

:j^ I stochastic background of gravitational waves has on pulsar timing parameters through 

C^ . the use of simulations and data from the millisecond pulsars PSR J0437-4715 and 

PSR J1713+0747. We show that the reported timing uncertainties become underesti- 
mated with increasing background amplitude by up to a factor of ~ 10 for a stochastic 
^ . gravitational-wave background amplitude of A = 5 x 10~^^, where A is the amplitude 

^^ ' of the characteristic strain spectrum at one- year gravitational wave periods. We find 

evidence for prominent low-frequency spectral leakage in simulated data sets includ- 
^^ ■ ing a stochastic gravitational-wave background. We use these simulations along with 

■^ ' independent Very Long Baseline Interferometry (VLBI) measurements of parallax to 

r — [ . set a 2-sigma upper limit oi A < 9.1 x 10"^*. We find that different supermassive 

f^ ■ black hole assembly scenarios do not have a significant effect on the calculated upper 

limits. We also test the effects that ultralow-frequency (10^^^-10~^ Hz) gravitational 
waves have on binary pulsar parameter measurements and find that the corruption of 
these parameters is less than those due to 10^^-10^^ Hz gravitational waves. 

Is^ 1 INTRODUCTION loss due t o gravitational radiation, as predicted by general 

rS . relativity (JTavlor fc WeisberdflQSa ). 

C^ ■ Pulsars are highly magnetized, rapidly rotating neutron Nearly three decades ago, it was first sh own that pul- 

stars that hav e spin periods rang i ng fro m seconds to mil- sars could be u sed to directly detect GWs (|Sazhinl Il978l : 

liseconds. See iLorimer &: Kramen i|2005l 1 for a full review iDetwe ilcr 1979) by using the pulsar timing residuals (fur- 

of pulsars and their applications. Pulsars can be separated ther di scussed in Section [21) to look for a specific GW sig- 

into two broad categories: normal pulsars and millisecond nature. iHellings fc Downs! (|l983l ') showed that by exploiting 

pulsars (MSPs). Normal pulsars are characterized by large the known correlations that a stochastic GW background 

periods and period derivatives (P ~ 1 s and P ~ 10"^'^, re- (GWB) would induce in the pulsar timing residuals, one 

spectively) and large magnetic fields {B ~ 10^^ G), whereas can place an upper limit on the GWB. This work introduces 

MSPs have smaller periods and period derivatives (P ~ 0.01 what is today known as the Hellings and Downs curve, which 

s and P ~ 10"^", respectively) and smaller magnetic fields is a main feature in many detection algorithms (see, e.g., 

(B ~ 10* G). Pulsars are extremely accurate clocks and, Ijenet et al.ll2005l : lAnholm et al.ll2009l 'l. The concept of a pul- 

specifically, MSPs exhibit much better timing stability as sar timing array (PTA) composed of the be st timed MSPs 

they do not show evidence of rotational instabilities such was first developed o ver two decades ago (|Romanil 1 19891 : 

as timing noise or glitches, which are prevalent in normal [Foster fc Backeilll99d ). Today there are three main PTAs 

pulsars. In fact, many MSPs rival atomic clock s in their in existence with the goal of GW detection using pulsars: 

fractional stability (see Figure 24 of lLorimerll2008l ) . Because the European Pulsar Timing Array (EPTA; Ijanssen et ajj 

of this and the fact that many MSPs are in binary systems, [2008), the North American Nan ohertz Observato ry for 

they are extensively used for tests of general relativity. The Gravitational waves (NANOGrav; Ijenet et all 2009 ), an d 



first evidence for the existence of gravitational waves (GWs) the Parkes Pulsar Timing Array (PPTA: lManchesteiil2008l V 



come from the precise timing of PSR B1913-I-16 in which the all of which are in collabora tion to form the I nternational 

decrease of the orbital period is in accordance with energy Pulsar Timing Array (IPTA: lHobbs et al.ll20ld ). 
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1.1 Pulsar timing precision and parameter 
uncertainties 

While the high timing precision of MSPs has allowed the 
precise measurement of many parameters and effects, pul- 
sar timing suffers from the inherent fact that the noise 
sources that contribute to the timing residuals are not well 
understood. To first order, our timing residuals are com- 
posed of noise that follows Gaussian statistics and that has 
a white power spectrum. For each time-of-arrival (TO A), 
the amount of receiver noise is estimated by the TOA un- 
certainty, as determined from a least-squares fit to a high 
signal-to-noise ratio template profile. It is, however, com- 
mon for these uncertainties to be underestimated or for ad- 
ditional sources of white noise to contribute to the timing 
residuals. One reason might be self-standarding: in the pro- 
cess of deriving an average TOA for an observation, the in- 
tegrated pulse profile is correlated with a template profile. 
If the template profile is generated from addition of many 
observations, then at some low level the noise from the ob- 
servation can correlate with the noise in the standard, lead- 
ing to an underestimation of the TOA uncertainty. Also, the 
cross-correlation method typically used has been shown to 
underestimate TO A uncertainties in the high-noise regime 
l|Hotan et al.ll2005l V Poor selection of standard profiles can 
underestimate the TOA uncertainty as well. All of these ef- 
fects cause the pulsar timing parameters to be consequently 
more uncertain than initially expected. Assuming these ef- 
fects scale with the uncertainties they affect, one can at- 
tempt to mitigate this underestimation by multiplication of 
the TOA uncertainties by the square root of the reduced x^ 
value of the timing residuals, as is commonly done in other 
fields. 

Underestimation of the white noise contribution is 
not the only challenge in assessing pulsar timing param- 
eters and their true uncertainties, though. The presence 
of non-white noise in pulsar timing residuals is very com- 
mon for young pulsars and, as timing baselines increase, 
starts to become important in timing of some MSP s as well 
llShannon fc Cordedl2010l : IVerbiest"et al. 2008; Spla ver et all 
[2OO5I). While the spectral properties and origins of non- white 
noise in MSPs are as yet not understood because the white 
noise is still too prominent to allow detailed analysis, it has 
already been shown that the least-squares fitting performed 
in pulsar timing reports underestimated uncertainties for 
the parameters in the timing fit and may even co rrupt the 
parameter values themselves (|Verbiest et al.ll2008r ). A tech- 
nique to improve the least-squares fitting process to mitigat e 
such effects has recently been proposed (jColes et al.ll201ll ). 

A third problem may be present in the form of signals 
that are part of a non-white spectrum - such as GWs. Not 
only does this type of signal affect the fitting process as 
described above, through introduction of power at a very 
specific frequency it has the potential to be fully absorbed 
in the timing signature of periodic parameters in the pulsar 
timing model, such as pulsar position, proper motion, par- 
allax or many binary parameters. On the one hand this can 
prove useful since independent measurement of parameters 
such as parallax (by, e.g., VLBI) can be combined with the 
pulsar timing measurement to place bounds on the am ount 
of corrupting noise present (see, e.g. iDeller et al.ir2008l ). but 
on the other hand it is worthwhile to assess the amount of 



influence such noise (especially predicted noise sources such 
as the GWB) can have on pulsar timing parameters. 



1.2 Pulsar timing and GWs 

Low frequency (10~^-10~^ Hz) GWs are expected from su- 
permassive black hole binary systems (SMBHBs), cosmic 
strings, and GWs from the big bang and inflationary era 
of the early universe. GWs from these sources can mani- 
fest themselves in different ways. Single nearby SMBHBs 
can pr oduce resolvable waves with periods on the order of 
years (|Wvithe fc Loebl [2003l ') . S MBHBs and cosmic strings 
can also produce GW bursts |Damour fc Vilenkinl I2OO1I : 
[Siemens et"aLl l2007l : iLeblond et al.ll2009l ) in which the du- 
ration of the GW signal is much less than the observation 
time. We also expect pulsar timing arrays to be sensitive to a 
stochastic background of unresolvable sources. This stochas- 
tic background can be described by a characteristic strain 
spectrum hc{,f) defined as a frequency-dependent power law 
(,Jenet et al.„200i l 



hc{f) = A 



f 



yr" 



(1) 



Here A is the dimensionless amplitude and a is the spectral 
index. The power in the GWB can then subsequently be 
written as 



Pif)-^js^^^(ff- 



(2) 



In accordance with many cosmological theories, one can also 
write r2gw(/), the energy density per unit logarithmic fre- 
quency interval, in terms of the characteristic strain spec- 
trum as 



"gw(/) 



2 Tv' 



fh4ff, 



(3) 



where Ho is the Hubble constant. This paper will focus 
primarily on a stochastic backgrou nd of SMBHB sources 
with a spectral index o f q = -2/3 JWvithe fc Loebl [2OO3I : 
[jaffe & Backer 2003; Enoki et al.i l2004l ). Recent work has 
shown that this simple power law does not hold true at 
higher frequencies (/ > 10~*) due to the discrete number 
of sources. At these frequenci es the characteristic strain is 
written as (|Sesana et al.ll2008l ) 



hcif) = ho 



1 + 



(4) 



where ho, fo, and 7 are model-dependent parameters based 
on specific SMBHB merger scenarios. The implications that 
this result will have on our simulations will be discussed in 
Section 14.31 where we test whether this broken power law 
will make a significant difference to methods that use peri- 
odic pulsar parameters with relatively high-frequency fitting 
functions. 

Although a direct detection of the GWB has not 
been made to date, methods using pulsar timing resid- 
uals have been used t o disprove proposed SMBHB sys- 
tems (| Jenet et al.1 120041') and to put limits on SMBHB co- 
alescence rates ( Wen et al.l 201 J) and on the dimensionless 
strain amplitude (IStinebring et aLJ Il990l : iJenet et ahl l2006l : 
Ivan Haasteren et al.l I2OIII ). The most s tringent publishe d 
upper limit to date is ^ < 1.1 x 10"^* l| Jenet et al.1 120061 ). 
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They use a statistic sensitive to a red spectrum. In this pa- 
per, we present a me thod which relies on the deriv ed timing 
IjVerbiest et al.ll200 9'l and VLBI dPeller et al.ll200 8l') parame- 
ters for individual pulsars. It involves adding simulated GWs 
to pulsar timing parameters and refitting for the pulsar tim- 
ing parameters to determine the effect of the induced GW. 
This idea of using VLBI and timing parameters to place 
limits on th e stochastic background was first introduced in 
IPeller et al.l ([2008). In that work, the kinematic distance ob- 
tained from timing measurements of Pb in the presence of a 
GWB is compared to the VLBI measurement of the distance 
to place an upper limit. 



1.3 This paper 

In this paper, we determine how strongly pulsar parame- 
ters can be corrupted by a GWB and red noise in gen- 
eral. We also use independent VLBI measurements to con- 
struct upper limits on the stochastic GWB. This is an 
important question because these corruptions shown here 
through simulations could be confirmed by more accurate 
VLBI measurements as we ll as other optical interfe rome- 
try missions such as GAIA JMignard fc Klioneij|2010l '). Fur- 
thermore, many relativistic tests in the strong field regime 
could be affected by using the corrupted values of the post- 
Keplerian parameters Pb and ijj. 

This paper is organized as follows: in Section [2] we de- 
scribe the basis of our simulations, including a brief review 
of GW effects on the pulsar timing residuals and fitting pro- 
cedures. In Section [3] we describe the observations used in 
our simulations. In Section |4] we present the effect of GWs 
on pulsar timing parameters and in Section [5] we summarise 
our findings. 



2 SIMULATION 

All simulations and timing analyses for this pa per were per- 
form ed using the Tempo2 software package (JHobbs et al.l 
120061 ). Full details of the timing m odel and algorithms em - 
ployed by Tempo2 can be found in lEdwards et al.l ()2006h . 



2.1 Pulsar timing theory 

Here we will only briefiy review the timing model and the 
standard fitting procedure. TOAs for a specific pulsar are 
produced at each epoch by a cross-correlation of the pulse 
profile with a low noise standard template. This procedure 
produces site-arrival-times (SATs), which are the times that 
the average pulse reached the telescope. The SATs are then 
converted to barycentric-arrival -times (BATs) using the So- 
lar System Ephemeris DE405 (|Standishl |2004| ) which gives 
accurate predictions of the locations and masses of Solar Sys- 
tem objects and the position of the telescope at any point in 
time. The timing model used in Tempo2 requires calculating 
the time of emission at the pulsar 

ir = if=-A0-Ais-AB, (5) 

where tf"" - A© = ff^^ is the BAT, A© includes all terms 
that go into converting the observed TOA to the solar 
system barycenter (SSB), (Roemer, Einstein, and Shapiro 
delays etc.), Ais includes delays due to passage through 



the interstellar medium, and Ab includes delays related 
to the pulsar's binary motion. Following the notation of 
I Edwards et al.l (|200y), the subscript "a" represents an ar- 
rival time and the "e" represents an emission time. In order 
to calculate the timing residuals, the pulse phase (/>(i) must 
be calculated and compared to the nearest integer value. 
The function </>(t) describes the time evolution of the pulse 
phase and is written in terms of a power series in time 



1 o 1 ^ 



(6) 



where At — t^^'^ — to is the difference in the time of emission 
from the pulsar and some reference epoch. From this it can 
be seen that fitting for the pulsar frequency v removes a 
linear offset from the timing residuals and fitting for the fre- 
quency derivative v removes a quadratic function. Although 
(l){t) is expressed as a power series in At, second order ap- 
proximations are very good for MSPs since they have neg- 
ligible i>. The fractional part of 4){t) is the "residual". The 
fitting procedure for this is a x^ minimisation where 



^ ,,2 Z^ 



(^(t») - N, 



(7) 



Here Ni is the nearest integer to 4>{ti), ai is the TOA uncer- 
tainty, and v is the pulse frequency. Since the fitting proce- 
dure is a x^ minimisation technique, the TEMP02-reported 
uncertainty is just the corresponding statistical uncertainty 
on th e reduced x^ St for a^ny given parameter ([Press et al.l 
Il992l ). 



2.2 Basic GW simulations 

Here we will give a brief review of the method through which 
Tempo2 adds a stochastic GWB to the pulsar timing resid- 
uals (see Hobbs et a,L (2009,) for full details). Tempo2 uses 
a time domain method to inject timing residuals caused by 
GWs. A stochastic background of GWs can be written as 
a sum of plane waves originating from random positions on 
the sky. The metric perturbation can then be written as 



hu 



Re 



JV-l 




- 


Y. ^A.-ie' 


(kj-X- 


-^jt) 


j=0 




. 



(8) 



and Uj are the position vector of the source, time, amplitude, 
wave number, and angular frequency of the jth GW source, 
respectively. It then follows that the residual induced by this 
sum of plane GWs is the fractional change in the frequency 
of the pulse rate integrated over time 
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Here repeated upper and lower indices indicate a sum over 
spatial coordinates, fcp is the unit vector in the direction of 
the pulsar, D is the distance from the pulsar to the Earth, 9j 
is the angle between the direction to the pulsar and the jth 
GW source, and Aimj is the transverse-traceless complex 
amplitude that takes the two GW polarisation states into 
account. 
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2.3 Extended GW simulations 

The basic algorithm used in this paper involves creat- 
ing 1,000 sets of induced timing residuals for each of 200 
GWB amplitudes for each pulsar, adding them to the ob- 
served residuals and refitting for the full timing model 
l|Verbiest et all |2008| . bOOgl ). Each amplitude is due to an 
ensemble average of 1,000 individual sources distributed 
isotropically on the sky. These sources are drawn from a 
distribution that follows Equation [T] The end result of this 
process is a distribution of the fitted parameters at 200 
GWB amplitudes distributed logarithmically in the range 
1 X 10~^^-1 X 10"^"^. This range was chosen because it in- 
cludes all previous upper limits and allows for lower back- 
grounds. These distributions were then analysed in several 
ways to calculate the effects of fitting for pulsar parame- 
ters when a GW signal is present in the timing residuals 
and to determine upper limits on the GWB amplitude. It is 
likely that there already is a GW signal from the stochastic 
background present in the pulsar timing data. Because we 
could be adding simulated GWs to those already existing in 
the data, any upper limits that we set could be, at worst, 
underestimated by a factor of \f2. 



the figures that the GWB can induce both high and low- 
frequency components in residuals. The problem with fitting 
for these parameters is that, in effect, one is fitting out a sig- 
nificant amount of the GW signal that is contributing to the 
residuals. It is important to note that although a stochas- 
tic GWB can induce large uncertainties in timing model 
parameters that are sinusoidally varying at yearly and half- 
yearly frequencies, a significant amount of the GW power 
absorbed comes from the lowest frequencies in the data due 
to spectral leakage. We have subsequently conducted addi- 
tional simulations in which we only include GWs in a very 
narrow frequency range around yearly and half-yearly fre- 
quencies. In these simulations we see that the corruption to 
the sinusoidally varying parameters is negligible even for a 
very large background amplitude of ^ = 5 x 10^^^. This 
leads to the conclusion that the absorbed power in the ear- 
lier simulations is indeed due to low-frequency spectral leak- 
age and not absorption of power at these higher frequencies. 
Nevertheless, absorption of high-frequency GW power still 
occurs, however, the steepness of the GWB power spectrum 
implies that for any realistic GWB this effect is fully covered 
up by the radiometer (white) noise. 



3 OBSERVATIONS 

We used data for the pulsars PSR J0437-4715 and PSR 
J1713-I-0747. These data were collected with the Parkes 64- 
m radio telescope at 20 cm wavelength. A full description 
of the data sets and observing systems can be found in 
IVerbiest et all (|2008l . [2009: 1 and references therein. These 
pulsars were chosen for three reasons. Firstly, the timing 
parallax measurements are consistent at the l-cr level with 
those from VLBI, giving us the ability to place upper lim- 
its on the background by comparing the timing and VLBI 
measurements. Secondly, both pulsars have low rms tim- 
ing residuals, which lead to the best upper limits that can 
be obtained using the methods described above. Finally, 
these data sets have a relatively long time span, allowing 
us to see effects of low-frequency GWs more readily in the 
timing residuals. A summary of the timing characteristics 
is given in Tables [l] and [2l For the VLBI parallax mea- 
surements in this study we use an LBA parallax for PSR 
J0437-4715 (iDeller et al.l l2008h and a V LBA parallax for 
PSR J1713+0747 (jChatteriee et al.ll2009l ). 



4 ANALYSIS 

4.1 Absorption of GW power in timing 
parameters 

Figure[T]illustrates how a GWB can be absorbed into certain 
periodic timing parameters for PSR J0437-4715. Figure [l(b)| 
shows the timing residuals with a GWB with amplitude of 
^ = 5 x 10~^^ added to the data. Figures 1(c) and 1(d) show 
the timing residuals using the GW- induced {A — ^x 10~^^) 
values of parallax and proper motion, and Pb, ^, and to, re- 
spectively. A very large GWB amplitude was used for illus- 
tration purposes, but this same effect takes place at smaller 
amplitudes and is just not as visible by eye. We can see from 



4.2 Impact on timing model parameters 

In this section we will construct confidence intervals for the 
fitted parameters at various GWB amplitudes and show that 
these confidence intervals become larger as the GWB ampli- 
tude increases. The method for calculating these confidence 
intervals is quite straightforward. The main simulation dis- 
cussed in Section 12.31 is used to obtain distributions of the 
given fitted parameters at each amplitude. As we have seen 
above, a GWB will cause excess low-frequency power to 
be absorbed into timing parameters. Examples of the sub- 
sequent corruption of the parameters can be seen in the his- 
tograms of Figure [2] Since this simulation injects GWs with 
random sky positions and random frequencies (within a cer- 
tain range and with a given spectrum), the subsequent cor- 
ruption of the fitted parameters caused by absorbing parts of 
the GW spectrum will result in a Gaussian distribution cen- 
tered around the unperturbed parameter values. It can be 
seen from Figure [2] that the FWHM of the distribution and 
therefore the confidence interval on the parameter increases 
with increasing GWB amplitude. This increased spread in 
potential parameter estimates is the GW-induced corrup- 
tion that we aim to quantify in this paper. It must be noted 
that at sufficiently high GWB amplitudes prominent low- 
frequency power will be visible in the timing residuals. To 
ensure phase-connection in our timing, we limited our simu- 
lations to a maximum GWB amplitude oflxlO"^"^. We have 
also used phase tracking to take the phase-wraps caused by 
the injected GWB into account in our simulations. We also 
note that towards the higher end of the simulated GWB 
amplitude range, sufficient levels of low-frequency noise will 
be present in the simulations to make the uncertainties re- 
tur ned from the standa rd least-squares fit unreliable (see 
e.g. IVerbiest et al.ll2008l ). However, since our analysis only 
uses the best-fit and not its uncertainty, this has no effect 
on our results. 

These histograms can be used to directly study the ef- 
fects that the presence of GWs has on the parameter esti- 
mates resulting from the fit. This is done by plotting the 
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ratio of the standard deviation of the fitted parameters and 
the TEMP02-reported unperturbed error on the parameters 
against GWB amplitude. The full list of fitted timing pa- 
rameters along with their values and uncertai nties are shown 
in Tables [1] and [2] for the unperturbed case (|Verbiest et alj 
|2008| . l2009i ) and the case where a GWB with amplitude at 
the current upper limit of A ~ 1 x 10^^'* was added to the 
data. The results are plotted in Figure [S] Firstly, it should 
be noted that the fits to spin frequency and frequency deriva- 
tive were by far the most affected by the GWB as can be 
seen in Tables [T] and [2] However, these parameters are not 
plotted here as this is expected since the spin frequency and 
frequency derivative fit out a linear trend and a quadratic 
trend respectively, thereby absorbing the lowest frequency 
(highest power) part of the GWB spectrum. We see that 
the orbital period derivative, Pb, is affected more for PSR 
J1713-H0747 than for PSR J0437-4715. This is likely due to 
the difference of orbital periods for these pulsars (5.7 and 
67.8 days for PSR J0437-4715 and PSR J1713-(-0747, re- 
spectively). The larger orbital period for PSR J1713+0747 
allows for lower frequency power to be absorbed resulting in 
a larger overall effect by Equation (2] 

Given that the parallax timing signature scales as 
D~^ cos P, where /3 is the ecliptic latitude and D is the pulsar 
distance, one would expect the effect of a GWB on the par- 
allax value will scale with the same factor. Using this scal- 
ing relation and the known ecliptic latitudes of the pulsars 
(30.87° for PSR J1713-h0747 and -67.87° for PSR J0437- 
4715), we expect that the effect of the GWB on the parallax 
should be a factor of ~ 3 times larger for PSR J0437-4715 
than for PSR J1713+0747. However, from Tables [l] and [2] 
we see that this ratio is 4.7. We also expect the effect of the 
GWB on the astrometric position to be stronger for PSR 
J0437-4715 than for J1713-I-0747 because of the lower eclip- 
tic of the former. However, we expect this difference to be 
less than that seen in our simulations. To find the origin of 
these discrepancies we created simulated data sets for the 
two pulsars with uniform sampling and bi-monthly observa- 
tions over 15 years. We then ran our GWB simulation on 
these data sets and found that the parallax corruption is 3 
times worse for PSR J0437-4715 than for PSR J1713-h0747 
and the astrometric position is slightly more affected for the 
former, as expected. The reason that the parallax and astro- 
metric position are not so heavily affected in the real data is 
because of the irregular sampling, which leads to correlations 
in the fitted parameters. To confirm this, we also simulated 
data sets with irregular spacing and compared their covari- 
ance matrices with those of the uniformly sampled data. We 
found that the correlations between parameters are lower for 
the evenly spaced data. Finally, it is also important to note 
that, for both pulsars, the increase in uncertainty of these 
parameters is approximately linear in GWB amplitude over 
the region of interest of 1 x 10"^^ < A<5x 10"". 



4.3 A limit on the GWB strength 

As discussed in Section [T] many papers have been published 
dealing with calculation of upper limits on the stochastic 
GWB amplitude. While some calculate this upper limit us- 
ing a method that is des i gned to detect the background 
(|van Haasteren et al.ll2009l . boiin . others use a statistic that 
is exclusively designed for putting limits on the background 



(|jenet et al.ll2006l ). For example, Jenet et al. use a statistic 
that is sensitive to a red power spectrum. By design, this 
method relies on the noise having a very white power spec- 
trum and thus any deviation from Gaussian noise is a prob- 
lem for this test. Since some non-Gaussian noise may occur 
in timing data, especially in long data sets, this method is 
limited in applicability. However, it has thus far produced 
the most stringent published limits. 

The method discussed in this paper is also not a suit- 
able candidate for detection but it does produce independent 
upper limits that are consistent with previously published 
limits. Before discussing specific methods. Figure |4] shows 
some general results of this method. To obtain these plots, 
we run the simulation as described in Section [2.31 and then 
take the output at each GWB amplitude and make a his- 
togram as shown in Figure [21 We calculate an upper limit 
using this scheme by finding the amplitude for which 95% of 
the simulations yield parameter values that lie outside of the 
formal 2-u errors on that parameter. Note here again that 
the proper motion is significantly affected for both pulsars. 
However, the orbital period derivative and astrometric po- 
sitions are strongly affected for PSR J1713+0747 and PSR 
J0437-4715 respectively, for reasons discussed in Section|4]2l 
While this analysis could be seen as providing a limit on the 
GWB amplitude, it depends on the timing value of these 
parameters, which is not fully independent of the simulated 
values. Therefore, we require an independent estimate of a 
timing parameter, as can be provided by VLBI, for example. 

We have shown that many astrometric parameters are 
affected by a GWB. However, when compa ring VLBI values 
(DcUcr ct al. 2008: 'Chat teriee et al.] |2009'') to timing values 
( Vcrbicst ct al. 2008 , 2009l 'l. it is readilv found that, with 
the exception of parallax, all VLBI parameter estimates are 
inconsistent with timing values at the many sigma level be- 
cause of calibrator uncertainty in the VLBI method. With 
the parallax method we repeat the procedure described in 
the previous section for the timing parallax, but instead of 
comparing against the unperturbed timing value and uncer- 
tainty, we use the VLBI parameter and parallax . To accom- 
plish this, we define an offset parameter 



A = 



TTVLBI — 7rT2| 
(7VLBI + CrT2 



(10) 



where ttvlbi and ■kt2 are the VLBI and Tempo2 simula- 
tion values of parallax, respectively and the denominator 
is the sum of the formal uncertainties on the parameters. 
The values of the VL BI parallaxes used in this paper are 
TV = 6.396 ±0.054 mas JPeller et al.ll2008l ) and n = 0.95i;;;!5^ 
mas JGhatteriee et al.l |2009') for PSR J0437-4715 and PSR 
J1713-I-0747, respectively. It is important to note that we 
compare the VLBI values with those published in the orig- 
inal timing models (IVerbiest et al.l 120081 . [2009|), where the 
reduced chi-squared was normalised through application of 
backend specific error-multiplication factors. While previous 
work has shown t hat a stochastic GWB will affect astro- 
metric par ameters (JGwinn et al.lll997l : lKopeikin et al.lll999l : 
|jaffell2004l ) measured by VLBI such as parallax, it is known 
that the angular deflections expected are on the order of the 
characteristic strain amplitude (see e.g. iBook fc FlanaganI 
I2OI0I ). Because typical amplitudes are ~ 10~^^ these effects 
are far less than the typical uncertainties on these astro- 
metric parameters. Thus, we can treat VLBI measurements 
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as unaffected by the stocfiastic GWB. To place a limit, we 
determine the GWB amplitude in which 95% of the reali- 
sations result in offset values in excess of 2. The factor of 2 
is to ensure a 5% false negative rate. The results from this 
method are shown in Table [S] It is clear from the table that 
the upper limits obtained from PSR J0437-4715 are more 
stringent than those obtained from PSR J1713-I-0747. This 
is to be expected, because the background does not cor- 
rupt the values of parallax as much for PSR J1713-I-0747, as 
shown in previous sections. Though using this method pro- 
duces upper limits that a re ~10 times less stringent than 
the best published limits (jjenet et al.l 120061 ) it provides an 
independent confirmation and is still constraining. 

4.4 Effects of the GWB spectral break 



Recent work (jSesana et al.l 120081 ) has shown that the stan- 
dard power law given in Equation [l] for the GW strain only 
holds for frequencies / < 10~^ Hz. Simulations were used to 
numerically formulate Equation |4] which describes the GW 
strain for frequencies / > 10~* Hz. GW detection strate- 
gies using pulsar timing are poised to detect the stochas- 
tic background at the lowest possible frequencies (~ 10~^ 
Hz), which is limited by the length of the data set, and 
thus still in the single power-law domain. The parameters 
/o,7, and /lo in Equation |4] depend on different scenarios 
for MB H formation and ev olution and can be found in Ta- 
ble 1 of lSesana et al.l (|2008h . F or clarity of our result s here, 
the four models used are VHM ("Volontcri ct al."200:j) , KBD 
ifKo ushiappas ct al. 2004), BVRhf (BcEclman ct al] |2006l ). 
and VHMhopk (|Lodato fc NataraianlbOOd ). 

A Tempo2 plugin was developed to test whether this 
broken power-law has any affect on pulsar timing parameters 
or our upper limit calculations described in the previous sec- 
tion. This plugin simulates a GWB using the simple power 
law spectrum for frequencies below the cutoff frequency of 
10~* Hz and then switch to the broken power law spectrum 
for frequencies greater than this cutoff. Since the param- 
eters are model dependent, the plugin allows the user to 
input which model to use thereby setting the values of /o 
and 7; however, /lo was simply determined by requiring that 
Equations [T] and 3] are continuous at the cutoff frequency. 
This step is required because producing an upper limit re- 
quires that we vary the value of the GW strain amplitude 
A and thereby ho. The strain spectrum for the above men- 
tioned MBH formation scenarios is plotted in Figure [S] It 
is clearly visible from the figure that there is a significant 
deviation from the power law strain spectrum. It also can be 
seen that the choice of model can make a considerable differ- 
ence in the strain spectrum. Figure [6] shows the strain of the 
individual sources used to create a stochastic background 
using the VHM model. Note the significant deviation from 
power law strain spectrum (dashed line). It should be noted 
here that we are using the average GW strain spectrum. In 
reality the spectral break comes from the small number of 
sources at higher frequencies, leading to a more jagged strain 
spectrum at these higher frequencies. While this method is 
sufficient for the work presented here, more realistic simula- 
tions will be very valuable for future studies of the stochastic 
GWB. However, such simulations are beyond the scope of 
this paper. The simulation was run on the PSR J0437-4715 
and PSR J1713-I-0747 data sets using the same method as 



described above for all four models. The results are sum- 
marised in Table [3] It is obvious that implementation of this 
broken power law does not significantly change the resulting 
upper limits. This is to be expected since the low-frequency 
GW power dominates the corruption of the timing parame- 
ters. In fact, the induced rms of a GWB at the frequencies 
corres ponding to the spectral b reak (~ 10~* Hz) is around 
10 ns (jSesana fc Vecchid 12010 '). which is a factor of ~ 20 
below our white noise level in both pulsars. We only include 
these results to conclusively show that the GWB spectral 
break has no effect on pulsar timing parameters. 



4.5 Effects of ultralow— frequency gravitational 
waves 



Previous work (|Kopeikinll 19971 ) has shown analytically that 
the secular variations of orbital period Pb and of the semima- 
jor axis a projected onto the line of sight (x = a sin t) , where 
L is the inclination angle, can be used to set upper limits on 
the energy density Ogw in the ultralow- frequency (10~^^- 
10~^ Hz) regime. This work derives upper limits in terms 
of the variance on the parameters Pb and x indicating that 
a stochastic GWB in the ultralow-frequency range would 
significantly corrupt these parameters. To see if this propo- 
sition is consistent with our simulations, we perform the 
same steps described in Section [4.11 except now we constrain 
the simulations to only the ultralow-frequency waves. Figure 
[7] show the results of these simulations. Although it is not 
shown in the figure, as in the normal frequency case there is 
no biasing of these parameter distributions. As shown in Fig- 
ure 7, the effects of ultra-low-frequency GWs is substantially 
smaller than that of low-frequency GWs. That notwith- 
standing, the data sets affected by ultra-low-frequency GWs 
were more severely affected by phase-wrapping, but this is 
primarily because the simulated frequency range (/ < 10~^) 
implies that to first order, the GWB inserts a strong linear 
slope into the timing residuals, which can easily cause phase 
wraps (i.e. timing residuals beyond one pulse period). The 
low-frequency background (/ > 10~^) however, introduces 
timing residuals that are more akin to a random walk and 
while these typically have a strong linear component as well, 
it is far less steep and therefore less likely to cause phase 
wrapping. The analysis on these two pulsars suggests that 
this method of placing upper limits on the GWB using the 
variances of Pb and x, respectively, is ineffective, or at best, 
far less constraining than simply using the variances in the 
low-frequency range. 



5 SUMMARY 

We have introduced a way to investigate the effects of low- 
frequency noise in the form of a stochastic background of 
GWs on pulsar timing parameters. While this work focuses 
particularly on the stochastic GWB, the methods for inves- 
tigating the noise/GW induced corruption of pulsar timing 
parameters could apply to any low- frequency noise source. 
This method involves a Monte-Carlo simulation over dif- 
ferent possible GWBs, adding them into the pulsar timing 
residuals, fitting for pulsar timing parameters, and produc- 
ing distributions of various GW-affected pulsar parameters. 
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These parameter distributions are then analysed to deter- 
mine the overall effect that the GWB has on pulsar timing. 
To summarise our results, we have: 



• Replicated the lVerbiest et alj (|2008l ) results of red noise 
leakage and how it affects parameter uncertainties, showing 
that some pulsar parameters could be underestimated by 
up to a factor of ~ 10 for a GWB with an amplitude of 

• Expanded on the above work with another pulsar (PSR 
J1713-I-0747), clarifying the importance of ecliptic latitude 
and the reliability of astrometric parameter measurements. 

• Shown that one can place upper limits on the stochastic 
GWB by comparing VLBI derived astrometric parameters 
to those obtained through pulsar timing. 

• Demonstrated practically that the spectral break (or 
its monochromatic components) are unlikely to affect pulsar 
timing before the SKA era. 

• Demonstrated that ultra-low frequency GWs have no 
obvious effect on pulsar timing parameters. 

• Demonstrated that even red noise at very low levels 
can leak into higher frequency terms and cause parameter 
corruptions. This lends a partial explanation for EFACs (a 
multiplicative error factor that is used to normalise the re- 
duced chi-squared in the fitting process) , even in a seemingly 
white data set. 

Future work could investigate these timing effects on a 
larger sample of pulsars, including normal (non-MSP) pul- 
sars to check for corruption of pulsar parameters though 
these effects should be smaller due to larger errors on the 
timing parameters. 
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Table 1. PSR J0437-4715 Timing Model Parameters 
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Parameter Name and Units Parameter Value Tempo2 Error" GWB Induced Error*" Error Ratii 

Fit and Data Set 

MJD range 50191.0-53819.2 

Number of TOAs 2847 

Observation length, T^^,^ (yr) 9.9 

rms timing residual (/xs) 0.199 

Measured Quantities 

Right ascension, a (J2000.0) 04 3715.8147635 3 57 19 

Declination, S (J2000.0) -47 15 08.624170 3 60 20 

Proper motion in a, fia (mas yr^^) 121.453 1 13 13 

Proper motion in S, /ig (mas yr~^) -71.457 1 13 13 

Annual parallax, n (mas) 6.65 7 98 14 

Dispersion measure, DM (cm"'^ pc) 2.6443 4 24 6 

Pulse frequency, u (ms) 173.68794618476804 3 2910 970 

Pulse frequency derivative, u (10"^^ s^) -1.7284079 3 585 195 

Orbital period, Pt (days) 5.74104646 108 216 2 

Orbital period derivative, Pb, (10^^^) 3.73 2 6 3 

Epoch of periastron passage, Tq (MJD) 52009.852429 582 582 1 

Projected semi-major axis, x (s) 3.36669708 11 11 1 

X (10~i*) 1.2'= 3= 27^= 9^= 

Longitude of periastron, a;o (deg) 1.2224 365 365 1 

Periastron advance, tj (deg yr-i) 0.01600 430 862 2 

Longitude of ascension, U (deg) 207.8 23 92 4 

Orbital inclination, i (deg) 137.58 6 24 4 

Companion mass, m2 (Mq) 0.25 1 11 

Orbital eccentricity, e (10-5) 1.9179 3 9 3 

Set Quantities 

Reference epoch for P, a, and S determination (MJD) 52005 

Reference epoch for DM determination (MJD) 53211 

" Given uncertainties are 1 a values in the last digits of the parameter values. 

^ 1 cr uncertainties for a simulated GWB with amplitude A = 1 X 10"^*. 

'^ Not part of original timing model in lVerbiest et al.l 1 120091 ). The full timing model including x was fitted in a separate simulation. 
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Table 2. PSR J1713+0747 Timing Model Parameters 



Parameter Name and Units Parameter Value Tempo2 Error" GWB Induced Error'' Error Rati( 

Fit and Data Set 

MJD range 49421.9-54546.8 

Number of TOAs 392 

Observation length, T^^,^ (yr) 14.0 

rms timing residual (/.ts) 0.198 

Measured Quantities 

Right ascension, a (J2000.0) 17 13 49.532628 1 6 6 

Declination, S ( J2000.0) +07 47 37.50165 3 12 4 

Proper motion in a, fia (mas yr~^) 4.924 5 35 7 

Proper motion in S, /ig (mas yr~^) -3.82 1 8 8 

Annual parallax, n (mas) 0.94 5 15 3 

Dispersion measure, DM (cm~'^ pc) 15.9915 1 3 3 

Pulse frequency u (Hz) 218.8118404414362 2 1786 893 

Pulse frequency derivative, u {W'^^ s^) -4.08379 2 134 67 

Orbital period, Pt (days) 67.825130963 9 27 3 

Orbital period derivative, Pb, (10"") 41 10 60 6 

Epoch of periastron passage, Tq (MJD) 54303.6328 4 8 2 

Projected semi-major axis, x (s) 32.3424236 2 4 2 

X (10~") -1.6^= 5"= 35° 7= 

Longitude of periastron, ujq (deg) 176.190 2 12 6 

Longitude of ascension, f! (deg) 67 9 63 7 

Orbital inclination, i (deg) 78.6 9 72 8 

Companion mass, rn2 (Mq) 0.20 2 12 6 

Orbital eccentricity, e (10-5) 7.4940 3 18 6 

Set Quantities 

Reference epoch for P, a, and 5 determination (MJD) 54312 

Reference epoch for DM determination (MJD) 54312 

" Given uncertainties are 1 a values in the last digits of the parameter values. 

' 1 cr uncertainties for a simulated GWB with amplitude A = 1 X 10"^*. 

'^ Not part of original timing model in lVerbiest et al.l 1 120091 ). The full timing model including x was fitted in a separate simulation. 



Table 3. Upper limits on the stochastic GWB amplitude (measured in units of 10 
models for pulsars PSR J0437-4715 and PSR J1713+0747. 



') at the 2-sigma level for different power law 



Spectrum PSR J0437-4715 PSR J1713-I-0747 

1.3 
1.3 
1.2 

1.2 
1.3 



a = -2/3 


0.91 


VHM 


0.91 


BVRhf 


0.91 


VHMhopk 


0.92 


KBD 


0.91 
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Figure 1. Plots of pulsar timing residuals of pulsar PSR J0437-4715 jVerbiest et al.ll200^ for different situations. Top left: no GWB 
added. We see here that the rnis is low and the noise is relatively white. Top right: post-fit residuals after adding a GWB with amplitude 
j4 = 5 X lO"'^^. Bottom left: no GWB added but using the GW-induced values of parallax and proper motion in our timing model. This 
shows strong periodic structure at periods ~ f yr. Bottom right: no GWB added but using GW-induced values of Pb, ui, and ui. 
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Figure 2. Distribution of the timing-derived parallax for PSR J0437-4715 as a result of fitting a timing model to 1000 realisations of 
residuals with a stochastic GWB added. Plotted here are the histograms (solid lines) of the fitted parameters as well as Gaussian curves 
(dashed lines) produced from the amplitude of the histogram, the mean of the distribution, and the standard deviation. It can be seen 
that a Gaussian fits the histogram very well. Here the histogram and Gaussian curve for the wider distribution is for a GWB amplitude 
of A ^ 5 X lO^"*^^ and the narrower distribution is for a GWB of A ~ 1 X 10^"'^^. 
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Figure 3. Plot of error ratio (ratio of standard deviation of distribution to the unperturbed error) vs. GWB amplitude for selected 
timing model parameters, excluding spin period and spindown (see Tables [T] and [2]l of PSR J0437-4715 and PSR J1713+0747. (Solid 
line: parallax, Dashed line: proper motion in right ascension. Dot-dashed line: proper motion in declination. Dotted line: right ascension. 
Dash-triple dot: declination.) It should be noted that the proper motion in right ascension and declination curves overlap for PSR 
J0437-4715. 
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Figure 4. Plot of detection significance vs. GWB amplit ude for PSRs PSR J0437-4715 and PSR J1713+0747 with parameters assumed 
to have the values listed in IVerbiest et al.l 1(2003, [2009J). The method for obtaining these plots is described in the text. These plots 
represent detection rates with false negative rates of 5% against GWB amplitude. Solid line: parallax, Dashed line: proper motion in 
right ascension, Dot-dashed line: proper motion in declination. Dotted line: right ascension. Dash-triple dot; declination.) 
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Figure 5. Plot of hc{f) against observed GW frequency. Here we show the naive power law spectrum from Equation [T] and the broken 
power spectrum from Equation |4] using different SMBH assembly models. (Thick dashed line: Normal power law spectrum, solid line: 
VHM model, dashed line: VHMhopk model, dash-dot line: KBD model, dotted line: BVRhf model.) 



< 
+ 



-9.5 







-8.5 



-7.5 



log,o[fg] (Hz) 



Figure 6. Plot of GW strain squared vs. GW frequency. The dotted line is the naive spectrum from Equation [T] for the GW strain 
spectrum. The points are individual SMBH binary sources emitting at random frequencies in the range 10~^-10~^ Hz generated by 
using Equation |4] and the VHM model. Note that there is a significant deviation in the spectrum around / > 10^** Hz. 
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Figure 7. Plot of error ratio vs. GWB amplitude for P^ and x in tlie normal (10^^—10^^ Hz) and ultralow (10^^^— 10^^ Hz) frequency 
ranges for PSR J0437-4715 and PSR J1713+0747. Tlic departure from a linear trend, seen at the highest simulated amplitudes, is caused 
by phase-wrapping that occurs when the simulated GW signature exceeds a few pulse periods. Such wrapping effectively randomises 
the results from the least-squares fit and produces outlier results in the Monte-Carlo iteration in which it occurs. This in turn distorts 
the statistics derived from that simulation and therefore corrupts the curves displayed here. We have used phase-tracking methods to 



avoid such wraps, but towards A ■ 
rendered useless. 
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